Lab-06: The Grammar of Maps

Which is your native?

Abstract

Part of the R for Artists and Designers course at the School of Foundation Studies, Srishti Manipal Institute of Art, Design, and Technology, Bangalore.

1 Goals

At the end of this Lab session, we should: - know the types and structures of spatial data and be able to work with them - understand the basics of modern spatial packages in R - be able to specify and download spatial data from the web, using R - plot static and interactive maps using ggplot, tmap and leaflet packages - add symbols and markers for places and regions of our own interest in these maps. - see directions for further work (e.g. maps + networks together)

2 Installation

rgdal is R’s interface to the “Geospatial Abstraction Library (GDAL)” which is used by other open source GIS packages such as QGIS and enables R to handle a broader range of spatial data formats.

2.1 Introduction to Maps in R

We will take small steps in making maps using just two of the several map making packages in R.

The steps we will use are:

1. Search for an area of interest ( E.g. using prettymapr or similar..)
2. Learn how to access spatial/map data using osmdata
3. Plot and dress up our map using osmplot, tmap and also with leaflet.
4. Create interactive maps with leaflet using a variety of map data providers. Note: tmap can also do interactive maps which we will explore also.

Bas. Onwards and Map-wards!!

2.2 God made me a BengaluR-kaR…I think

Let’s get BLR data into R and see if we can plot an area of interest. Then we can order on Swiggy…never mind.

Where is my home? Specify a “bounding box” first, using a rough longitude latitude info, or using a place name and searching for the long/lat:

# BLR Bounding Box
bbox <- osmplotr::get_bbox(c(77.56,12.93,77.63,12.96))
bbox_l <- osmdata::getbb("Bangalore, India")
bbox_p <- prettymapr::searchbbox("Bangalore")

## Using default API key for pickpoint.io. If batch geocoding, please get your own (free) API key at https://app.pickpoint.io/sign-up

bbox

##     min   max
## x 77.56 77.63
## y 12.93 12.96

bbox_l

##        min      max
## x 77.46010 77.78405
## y 12.83401 13.14366

bbox_p # identical with bbox_l

##        min      max
## x 77.46010 77.78405
## y 12.83401 13.14366

# Get Map data

dat_B <- extract_osm_objects (key = "building", bbox = bbox) 

## Issuing query to Overpass API ...

## Rate limit: 0

## Query complete!

## converting OSM data to sf format

dat_H <- extract_osm_objects (key = 'highway', bbox = bbox)

## Issuing query to Overpass API ...

## Rate limit: 0

## Query complete!

## converting OSM data to sf format

dat_P <- extract_osm_objects (key = 'park', bbox = bbox)

## Issuing query to Overpass API ...

## Rate limit: 0

## Query complete!

## converting OSM data to sf format

dat_G <- extract_osm_objects (key = 'landuse', value = 'grass', bbox = bbox)

## Issuing query to Overpass API ...

## Rate limit: 0

## Query complete!

## converting OSM data to sf format

dat_T <- extract_osm_objects (key = 'natural', value = 'tree', bbox = bbox)

## Issuing query to Overpass API ...

## Rate limit: 0

## Query complete!

## converting OSM data to sf format

#Useful keys include building, waterway, natural, grass, park, amenity, shop, boundary, and highway.

2.3 My first Map in R

blr_map <- osm_basemap(bbox = bbox, bg = "gray20") %>%
  add_osm_objects(., dat_B, col = "gray40") %>% 
  add_osm_objects(., dat_H, col = "gray80") %>% 
  add_osm_objects(., dat_G, col = "darkseagreen1") %>% 
  add_osm_objects(., dat_P, col = "darkseagreen") %>% 
  add_osm_objects(., dat_T, col = "green")

print_osm_map(blr_map)
#print_osm_map (map, filename = "Blr.jpeg", device = "jpeg", width = 10, units = "cm")
#print_osm_map (map, filename = "map.png", width = 10, units = "in", dpi = 300)

blr_map_2 <- osm_basemap(bbox = bbox, bg = "gray20") %>%
  add_osm_objects(., dat_B, col = "gray40") %>% 
  add_osm_objects(., dat_H, col = "gray80") %>% 
  add_osm_objects(., dat_G, col = "darkseagreen1") %>% 
  add_osm_objects(., dat_P, col = "darkseagreen") %>% 
  add_osm_objects(., dat_T, col = "green")

print_osm_map(blr_map)

class(dat_B)

## [1] "sf"         "data.frame"

names(dat_B)

##  [1] "osm_id"                      "name"                       
##  [3] "Friary"                      "addr.city"                  
##  [5] "addr.district"               "addr.full"                  
##  [7] "addr.housename"              "addr.housenumber"           
##  [9] "addr.postcode"               "addr.state"                 
## [11] "addr.street"                 "addr.suburb"                
## [13] "air_conditioning"            "alt_name"                   
## [15] "amenity"                     "area"                       
## [17] "atm"                         "barrier"                    
## [19] "brand"                       "brand.wikidata"             
## [21] "brand.wikipedia"             "building"                   
## [23] "building.colour"             "building.levels"            
## [25] "building.levels.underground" "bus"                        
## [27] "capacity"                    "clothes"                    
## [29] "club"                        "construction"               
## [31] "cuisine"                     "delivery"                   
## [33] "denomination"                "description"                
## [35] "designation"                 "ele"                        
## [37] "email"                       "emergency"                  
## [39] "government"                  "healthcare"                 
## [41] "healthcare.speciality"       "height"                     
## [43] "internet_access"             "landmark"                   
## [45] "landuse"                     "layer"                      
## [47] "leisure"                     "level"                      
## [49] "lit"                         "man_made"                   
## [51] "max_age"                     "min_age"                    
## [53] "name.en"                     "name.kn"                    
## [55] "note"                        "office"                     
## [57] "official_name"               "opening_hours"              
## [59] "operator"                    "operator.type"              
## [61] "organic"                     "outdoor_seating"            
## [63] "payment.cash"                "payment.debit_cards"        
## [65] "phone"                       "public_transport"           
## [67] "religion"                    "screen"                     
## [69] "service"                     "service_times"              
## [71] "shop"                        "smoking"                    
## [73] "source"                      "takeaway"                   
## [75] "tourism"                     "website"                    
## [77] "wheelchair"                  "wikidata"                   
## [79] "wikipedia"                   "geometry"

class(dat_B$geometry)

## [1] "sfc_POLYGON" "sfc"

nrow(dat_B)

## [1] 39699

2.4 Adding my “area” to the map

We can create areas of interest around the map. Say based on metro stations or your restaurants etc.

my_area <- 
  bbox %>% 
  zoombbox(factor = 16, offset = c(0,0)) 
my_area

##        min      max
## x 77.59281 77.59719
## y 12.94406 12.94594

bbox

##     min   max
## x 77.56 77.63
## y 12.93 12.96

my_area_in_blr <- 
  sp::SpatialPoints(coords = cbind(
  c(my_area["x", "min"], my_area["x", "max"],
    my_area["x", "max"], my_area["x", "min"]),
  c(my_area["y", "min"], my_area["y", "min"],
    my_area["y", "max"], my_area["y", "max"])
        )
  )
my_area_in_blr

## class       : SpatialPoints 
## features    : 4 
## extent      : 77.59281, 77.59719, 12.94406, 12.94594  (xmin, xmax, ymin, ymax)
## crs         : NA

blr_map %>% add_osm_groups(
  .,
  dat_B,
  groups = my_area_in_blr,
  cols = "yellow",
  bg = "gray40",
  colmat = FALSE
) %>%
  add_osm_objects(., dat_B, col = "red", size = 0.2) %>% 
  print_osm_map(.)

2.4.0.0.1 END OF CLASS 1

2.5 Using rnaturalearth and tmap

data("World")
data("metro")

head(metro, n = 3)

## Simple feature collection with 3 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 3.04197 ymin: -8.83682 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS:  WGS 84
##       name             name_long iso_a3 pop1950 pop1960 pop1970 pop1980 pop1990
## 2    Kabul                 Kabul    AFG  170784  285352  471891  977824 1549320
## 8  Algiers El Djazair  (Algiers)    DZA  516450  871636 1281127 1621442 1797068
## 13  Luanda                Luanda    AGO  138413  219427  459225  771349 1390240
##    pop2000 pop2010 pop2020  pop2030                  geometry
## 2  2401109 3722320 5721697  8279607 POINT (69.17246 34.52889)
## 8  2140577 2432023 2835218  3404575   POINT (3.04197 36.7525)
## 13 2591388 4508434 6836849 10428756 POINT (13.23432 -8.83682)

tmap_mode("plot")

## tmap mode set to plotting

tm_shape(World) +
    tm_polygons("HPI") +
  
tm_shape(metro) + 
  tm_bubbles(size = "pop2020", col = "red")

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(World) +
    tm_polygons("HPI") +
tm_shape(metro) + 
  tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")

## Legend for symbol sizes not available in view mode.

india <- 
  ne_states(country =  "india", 
            geounit = "india", 
            returnclass = "sf")

india_neighbours <- 
  ne_states(country= (c("india", "sri lanka", "pakistan", "afghanistan", "nepal","bangladesh")))

names(india)

##  [1] "featurecla" "scalerank"  "adm1_code"  "diss_me"    "iso_3166_2"
##  [6] "wikipedia"  "iso_a2"     "adm0_sr"    "name"       "name_alt"  
## [11] "name_local" "type"       "type_en"    "code_local" "code_hasc" 
## [16] "note"       "hasc_maybe" "region"     "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me"   "datarank"   "abbrev"     "postal"    
## [26] "area_sqkm"  "sameascity" "labelrank"  "name_len"   "mapcolor9" 
## [31] "mapcolor13" "fips"       "fips_alt"   "woe_id"     "woe_label" 
## [36] "woe_name"   "latitude"   "longitude"  "sov_a3"     "adm0_a3"   
## [41] "adm0_label" "admin"      "geonunit"   "gu_a3"      "gn_id"     
## [46] "gn_name"    "gns_id"     "gns_name"   "gn_level"   "gn_region" 
## [51] "gn_a1_code" "region_sub" "sub_code"   "gns_level"  "gns_lang"  
## [56] "gns_adm1"   "gns_region" "min_label"  "max_label"  "min_zoom"  
## [61] "wikidataid" "name_ar"    "name_bn"    "name_de"    "name_en"   
## [66] "name_es"    "name_fr"    "name_el"    "name_hi"    "name_hu"   
## [71] "name_id"    "name_it"    "name_ja"    "name_ko"    "name_nl"   
## [76] "name_pl"    "name_pt"    "name_ru"    "name_sv"    "name_tr"   
## [81] "name_vi"    "name_zh"    "ne_id"      "geometry"

names(india_neighbours)

##  [1] "featurecla" "scalerank"  "adm1_code"  "diss_me"    "iso_3166_2"
##  [6] "wikipedia"  "iso_a2"     "adm0_sr"    "name"       "name_alt"  
## [11] "name_local" "type"       "type_en"    "code_local" "code_hasc" 
## [16] "note"       "hasc_maybe" "region"     "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me"   "datarank"   "abbrev"     "postal"    
## [26] "area_sqkm"  "sameascity" "labelrank"  "name_len"   "mapcolor9" 
## [31] "mapcolor13" "fips"       "fips_alt"   "woe_id"     "woe_label" 
## [36] "woe_name"   "latitude"   "longitude"  "sov_a3"     "adm0_a3"   
## [41] "adm0_label" "admin"      "geonunit"   "gu_a3"      "gn_id"     
## [46] "gn_name"    "gns_id"     "gns_name"   "gn_level"   "gn_region" 
## [51] "gn_a1_code" "region_sub" "sub_code"   "gns_level"  "gns_lang"  
## [56] "gns_adm1"   "gns_region" "min_label"  "max_label"  "min_zoom"  
## [61] "wikidataid" "name_ar"    "name_bn"    "name_de"    "name_en"   
## [66] "name_es"    "name_fr"    "name_el"    "name_hi"    "name_hu"   
## [71] "name_id"    "name_it"    "name_ja"    "name_ko"    "name_nl"   
## [76] "name_pl"    "name_pt"    "name_ru"    "name_sv"    "name_tr"   
## [81] "name_vi"    "name_zh"    "ne_id"

2.6 Subsetting Spatial data by attributes

#ind <- metro$iso_a3 == "IND"
#metro_ind1 <- metro[ind,]
metro_ind2 <- subset(metro, iso_a3 == "IND")
metro_neighbours <- metro %>% dplyr::filter(iso_a3 %in% c("IND","PAK", "LKA", "BGD","NPL"))

2.7 Map 1

tm_shape(World %>% dplyr::filter(iso_a3 %in% c("IND", "AFG", "PAK", "NPL", "BGD", "LKA"))) + 
  tm_borders() + 
tm_shape(india) + 
  tm_polygons("name") + 
  tm_shape(metro_ind2)+
  tm_dots(size = "pop2020") + 
  tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
  tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "right") + 
  tm_compass(position = c("right", "top")) + 
  tm_scale_bar(position = "left") + 
  tmap_style("watercolor") #cobalt #gray #white #col_blind #beaver #classic #watercolor #albatross #bw

## tmap style set to "watercolor"

## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic"

## Credits not supported in view mode.

## Compass not supported in view mode.

## Warning: Number of levels of the variable "name" is 35, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 35) in the layer function to show all levels.

## Legend for symbol sizes not available in view mode.

2.8 Map 2

tm_shape(india_neighbours) + 
  tm_polygons("name") + 
  tm_shape(metro_neighbours) +
  tm_dots(size = "pop2020") + 
  tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
  tmap_options(max.categories = 10) + 
  tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "center")

## Credits not supported in view mode.

## Warning: Number of levels of the variable "name" is 112, which is
## larger than max.categories (which is 10), so levels are combined. Set
## tmap_options(max.categories = 112) in the layer function to show all levels.

## Legend for symbol sizes not available in view mode.

tmap_mode("plot")

## tmap mode set to plotting

tm_basemap("Stamen.Watercolor") +
  tm_shape(metro, bbox = "India") + 
  tm_dots(col = "red", 
          # user-chosen group name for layers
          group = "Metropolitan Areas") +
  tm_shape(World) + 
  tm_borders() + 
  tm_tiles(server = "Stamen.TonerLabels",group = "Labels") + # ADDS LABELS!!!
  tm_graticules()

2.9 Scope and Packages for Exploration!!

2.9.1 sfnetworks

2.9.2 mapsf

2.9.3 ggspatial

2.10 Resources

1. Emine Fidan, Guide to Creating Interactive Maps in R